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We present an improved lattice Boltzmann model for high-speed compressible flows. The 
model is composed of a discrete-velocity model by Kataoka and Tsutahara [Phys. Rev. 
E 69, 056702 (2004)] and an appropriate finite-difference scheme combined with an ad- 
ditional dissipation term. With the dissipation term parameters in the model can be 
flexibly chosen so that the von Neumann stability condition is satisfied. The influence of 
the various model parameters on the numerical stability is analyzed and some reference 
values of parameter are suggested. The new scheme works for both subsonic and super- 
sonic flows with a Mach number up to 30 (or higher), which is validated by well-known 
benchmark tests. Simulations on Riemann problems with very high ratios (1000 : 1) 
of pressure and density also show good accuracy and stability. Successful recovering of 
regular and double Mach shock reflections shows the potential application of the lattice 
Boltzmann model to fluid systems where non-equilibrium processes are intrinsic. The 
new scheme for stability can be easily extended to other lattice Boltzmann models. 

Keywords: Lattice Boltzmann; high-speed compressible flow; von Neumann Analysis; 
shock. 



1. Introduction 

High-speed compressible flow with shocks plays an important role in various fields, 
such as explosion physics, aeronautics, etc. EfRcient simulation of such a system 
is interesting and challenging. The traditional method is based on a set of macro- 
scopic Euler equations resolved by the Finite-element or Finite-volume schemes, 
wher e the artificial viscosity is applied or the Riemann solver is used to capture the 
shociffEEl. According to the gas kinetic theory, a set of Euler equations describes a 
system being at equilibrium. For a system with shocks, the non-equilibrium behavior 
is intrinsic, so a scheme based on the fundamental kinetic theory is to be preferred. 
As a new approach to fluid dynamics, the lattice Boltzmann (LB) methoc3 solves 
the fully discrete Boltzmann equation by using an appropriate difference scheme to 
the temporal and spatial derivatives of the distribution function /^(x, t), where x 
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and t are the position and time, respectively, and the index i corresponds to the 
i-th discrete velocity. It recovers the desired macroscopic equations in the hydrody- 
namic limit and has the potential to fill the gap between continuum description and 
molecular dynamics'' . Besides the traditional LB originating from the lattice gas cel- 
lular automata^ 7 8 9 10^ other versions such as finite-difference(FD)-'^-'^ ^'^ 14 15]^ 
finite- volume(FVjJ^, and finite-element (Fe|I^, etc have also been developed un- 
der the same framework. Among these works, developing LB mo dels fo r high-speed 
compressible flows has long been attempted by different 

authorA^llSI. 

Among the 

existing models for two-dimesnional compressible fluids, the one by Kataoka and 
Tsutahara(KT|lSl]^a,s a simple and rigorous theoretical background. It takes flexible 
ratio of specific-heat and is superior in computational efficiency because the t otal 
number of its discrete velocity is reduced to 9. But similar to previous LB modelJ^, 
the numerical stability problem remains one of the few blocks for its practical sim- 
ulation to high-Mach-number compressible flows. In this paper we present a new 
scheme based on the original discrete- velocity-model (DVM) by KT and an appro- 
priate finite-difference scheme combined an additional dissipation term. With the 
new scheme fluid systems with high-Mach-number and/or high ratios of pressure 
and density can be successfully simulated. 

This paper is organized as follows. In section 2 the original discrete-velocity- 
model by KT is briefly reviewed and an alternative FD scheme is proposed for later 
analysis and using. A von Neumann stability analysis is performed in section 3, from 
which solutions to improve the numerical stability can be found. Several benchmark 
tests are used to validate the proposed scheme in section 4. Section 5 concludes the 
present paper. 



2. Description of the DVM and FD scheme 

The LB equation with the Bhatanger-Gross-Krook approximatiorPQl reads, 

dfi dfi 1 g„ 

where f^'^ is the discrete version of the local equilibrium distribution function; r the 
relaxation time; index a = 1, 2, 3 corresponding to x, y, and z, respectively; and 
Vi the i-th discrete velocity, i — 0, A^ — 1; A^is the total number of the discrete 
velocity. Under the hydrodynamic limit the LB equation is required to describe the 
following Euler equations, 

dp ^ d{pua) ^ ^ 



dt dXa 
djpuq) _^ d{puaup) ^ dP 
dt dxp dxa 

dp{bRT + ul) ^ dpUq{bRT + ul) + 2Puq 
dt dx/3 



= 0, (2) 

= 0, 
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where p, u, T, P (= pRT) are the hydrodynamic density, flow velocity, temperature 
and pressure, respectively, and R is the specific gas constant, h relates to the specific- 
heat ratio 7 as follows, b — 2/(7 — 1). The following constraints are imposed on the 
moments of fl'^ and /i, 

N-l N-1 

i=0 i=0 

N-l N-l 

pUa = ^ fi'^Via = ^ fiVia, (4) 

1=0 i=0 

N-l N-l 

p{bRT + ul)^Yl ft'^o. + ^l) = E + (5) 

i=0 1=0 

N-l 

PSaf3 + pUaUp = E ffviaViji, (6) 
i=0 



N-l 

p[{b + 2)RT + = E + )^^», (7) 

i=0 



where rji is another variable introduced to make specific-heat ratio flexibl^. 

Equation ([1]) may be written in non-dimensional form by using a characteristic 
flow length scale L, reference speed Cr and density Pr- Two reference time scales 
are used, tc to represent the time between particle collisions and L/e^ to present a 
characteristic flow time. The resulting non-dimensional equation is 

where the caret symbols are used to denote non-dimensional quantities ijia — Via/cr, 
t = tCr/L , f — r/tc, and fi — fi/pr- The parameter e — tcCr/L is the Knudsen 
number which may be interpreted as either the ratio of collision time to flow time 
or as the ratio of mean free path to the characteristic flow length. We will not use 
the caret notation further but will assume that the equation are in non-dimensional 
form henceforth. 

In the two-dimensional case, the KT discrete velocity model has nine compo- 
nents. It reads 

r (o,o),z = o 

{v^l,v,2)^< ci(cos(2(^),sin(2(^)),z = l,2,3,4 (9) 
[ C2 (cos7r(i±i + i),cos7r(^ + 3)) = 5,6,7,8 

^In a practical system, the ratio 7 provides information on the internal degrees of freedom of 
molecules. For example, 7 has a certain well-known value for an ideal, monatomic gas (like helium), 
and is different for diatomic molecules like those that make up most of the atmosphere. To formulate 
the DVM, the discretization and contribution of the internal degrees of freedoms of the molecules 
are represented by the constraints ((Sjl and l[7)l. 
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Fig. 1. (Left above) Schematic figure of the discrete velocity model. 

Fig. 2. (Right above) Effects of the dissipation term. Parameters used are p = 1.0, T = 1.0, 
ui = 10.0, U2 = 0.0, the remaining constants are set as ci = 10, C2 = 20, tiq = 10, Ai = aAx/c2 



Vi = 



i 





1,2,. 



(10) 



A schematic figure of the distribution of the discrete velocities is shown in Fig.l, 
where ci and C2 are constants which should not depart faraway from the flow velocity 
u and C2 is generally chosen 1.0 ^ 3.0 times of ci. 

The local equilibrium distribution function is computed by 



= p{Ai + BiViaUa + DiUaViaUi3Vii3), 1 ^ 0, 1, 2, , 



where 



A, = 



^-^T, i = 

1 



4(cf- 



i{ci-ci) 



{b-2 




i = 1,2,3,4 
i = 5,6,7,8 



(11) 



(12) 



B,, 



0, i = 

-4 + {b+2)T+ul 

-cf + {b+2)T+ul 
■^ciici-ci) ■ 



1,2,3,4 
5,6,7,8 



D, 




(13) 



It is clear that rjo, ci and C2 are independent parameters in this DVM and the 
value of ?7o influences the discrete local equilibrium distribution function f^"^ via 
the expansion coefficient Ai. The combination of the above DVM and the general 
FD scheme with first-order forward in time and second-order upwinding in space 
composes the original FDLB model by KT. The FDLB by KT has been validated 
via the Riemann problem in subsonic fiowJ^^. In a LB simulation the discretization 
in time and space introduces unphysical waves, and the collision term introduces 
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a physical dissipation when the system deviates from the local equilibrium. If the 
physical dissipation is strong enough so that the unphysical oscillations are not to 
be amplified in the simulation procedure, we will have no instability problem. The 
original LB model by KT is not stable when the Mach number M exceeds 1^^, which 
shows that an additional dissipation term is needed in such cases. To make practical 
the LB simulation to the supersonic flows, we propose an alternative FD scheme in 
the following part of this section. The proposed FD scheme will be combined with 
an additional dissipation term to overcome the numerical instability problem in the 
next section. 

We use the usual first-order forward scheme in time. Since all the quantities are 
now non-dimensional, to simplify the following analysis, the time step At is set to 
be numerically equal to the Knudsen number e. Thus, from Eq.® we have 

/,(x, i + At) - /,(x, t) + v^^^I^^At = i [/;«(x, t) - /,(x, t)] . (14) 

In Ea.(fH|) the spatial derivative dfi/dx can be calculated by 

^„ dh PMx + Ax, t) + (1 - 2/?)/,(x, t) - (1 - P)nix - Ax, t) 
If ... > 0, — = ; (15) 

„ (1 - P)Mx + Ax, t) ~ (1 - 2p)n{x, t) ~ PU{x - Ax, t) 

If < 0, — = . (16) 

ox Ax 



In Eqs. fTS)) and (|T6|) . < /? < 0.5. If /? takes zero, then they are not other than the 
first order upwind scheme in space; if /? takes 0.5, they recover to the general central 
difference scheme, dfi/dy can be calculated in a similar way. Actually, Eas. (fT5| and 
(jl6p can be rewritten as 

„ f,{x,t) - f,{x - Ax,t) 

, l3Ax[f,{x + Ax, t) + Mx- Ax, t) - 2Mx, t)] 



A: 



2 



lf..<0, 9f,^ Mx + Ax,t)-Mx,t) ^^^^ 
PAxlMx + Ax, t) + Mx- Ax, t) - 2Mx, t)] 



A; 



2 



The second terms in the right-hand-side of Eqs.(fT7| and (fT8|) can be regarded as 
some kind of artificial viscosities which are used to reduce some unphysical phe- 
nomena such as wall-heatinj^, but they are not enough to be effectively improve 
the stability of LB simulation, which means additional dissipation term is needed 
for a practical LB simulation. In the following sections the parameter (3 is chosen 
to be 0.25 if not particularly stated. 
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3. von Neumann Analysis 

The stab ility problem of LB has been addressed and attempted for som e years 
| 4 | 19 | 21 | 23 24 25 26 27 Among them, the the entropic LB met hocP^ l^ tries to 
make the scheme to follow the iJ-theorem; The FIX-UP methocPSEi] based on 
the standard BGK scheme, uses a third order equilibrium distribution function and 
a self-adapting updating parameter to ayoid negatiyeness of the mass distribution 
function. Flux limiter techniques are used to enhance the stability of FDLB by 
Sofonea, et aPSl. Adding minimal dissipation locally to improve stability is also 
suggested by Brownlee, et a 1221, but there such an approach is not explicitly dis- 
cussed. All the above mentioned attempts are for low Mach number flows. In this 
paper we focus mainly on h igh speed flows. 

Following Seta, et aP, in this paper we resort to the von Neumann stability 
analysis to compose a stable LB scheme where the additional dissipation is effective 
and minimal. The following analysis is based on the FD scheme shown in Eqs. pT]) 
and (|18p . In the von Neumann analysis the solution of finite-difference equation is 
written as the familiar Fourier series, and the numerical stability is evaluated by the 
magnitude of eigenvalues of an amplification matrix. The small perturbation A/j is 
defined as /i(x, i) = A/i(x, t) + /°, where /° is the global equilibrium distribution 
function and is a constant which does not vary in space or time and depends only 
on the mean density, velocity and temperature. From Eq. (|14p we can obtain 



A/,(x, t + Ai) - A/,(x, i) -I- «,„?^At = ^ 



(19) 



The perturbation part Afi{'x.,t) may be written as series of complex exponents, 
A/i(x, t) — i^|exp(ik • x), where is an amplitude at grid point x and time t, i 
is an imaginary unit, and is the wave number of sine wave in the domain with 
the highest resolution I/Axq. From Eqs. p9|) we obtain FI+^^ = G^^jFj , where 
Gij is a matrix being used to assess amplification rate of F* per time step At. If 
the maximum of the eigenvalues of the amplification matrix satisfies the condition, 
max|a;| < 1, for all wave numbers, the FD scheme is surely stable, where lu is 
the eigenvalue of the amplification matrix. This is the von Neumann condition for 
stability. 

The amplification matrix G can be written as following. 



where 



Pexp{ikaAxa) + (1 — 2/3) — (1 — /3)exp(— ifccAxa), if Via > 0; 
(1 — P)e'}q){ikaAxa) — (1 — 2/3) — /3cxp(— ifcaAxo,), if Via < 0. 



(21) 



Se veral rese archers have analyzed the stability of the incompressible LB 
model JTTO^Ql ^ 

it is found that there is not a single wave-number being always 
the most unstable. For the 2D DVM by KT G is a matrix with 9x9 elements. 
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Every element is related to the macroscopical variables (density, temperature, ve- 
locities), discrete velocities and other constants, so it is difficult to analyze with 
explicit expressions. We resort to the software, Mathematica-5. 

In order to simulate high-speed flows, we introduce the following dissipation 
term to the LB equation, 

/,(x,i+At)-/,(x,t)+z;„^^^At-A, ^ ^^|^Ai = ^[/^''(x, t)-/,(x, t)](22) 

where Xi is a small number not varying in space or time. The second-order derivative 
^ al^'*^ can be calculated by the central difference scheme. Then Gij can be written 
as 

dfP v^c.^t ' 2-^2cos(fc^x^ 

Obviously, in Eq. (j23p the last term is required to improve the numerical stability. 
How to chose the A; is the key problem here. It will not be effective if too small and 
will result in too additional errors if too large. To get some indication we look back 
to the last terms in Eqs. fTT)) and which are regarded as artificial viscosities to 
reduce the numerical wall-heating phenomena. To simplify the discussion, we choose 
always Aa; = Ay. Indicative analysis and numerical tests show that we can choose 
Xi around the following way, 




Xi = <^ ci Ax/lO, i = 1, 2, 3, 4 . (24) 



Now we show some results of von Neumann analysis by Mathematica-5 to get 
a more complete understanding of the stability condition. We will show only the 
results for high-Mach-number flows where the instability problem is generally much 
more pronounced and previous LB models fail to work. The results will be shown by 
figures with curves for the maximum eigenvalue jwlmaa: of G versus kAx. The wave 
number k is discretized from to tt with 30 steps. Figure 2 shows a comparison 
between the two cases, with and without the additional dissipation term, where the 
macroscopic variables are chosen as p = 1.0, T = 1.0, ui = 10.0, M2 = 0.0, and the 
constants in Eqs.® and (fTO|) are set as ci — 10, C2 — 20, rjo = 10. Coefficient a 
in the inset of the figure is a new constant introduced to control the time step in 
the following way. At = a Ax/ €2- For the two sets of results shown in the figure, 
it is clear that the dissipation term can significantly decrease the the maximum 
eigenvalue \uj\max from being larger than to be smaller than 1 for appropriately 
given time step. 

It is interesting to investigate the effects of various parameters (physical quan- 
tities and model constants in Eqs.Q and (fTO|) ) on the numerical stability. Fig. 3 
shows the a comparison of two cases: the first one is /3 = 0.25 with switching on the 
additional dissipation and the second is /? = with switching off of the additional 
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N, 




0.0 0.5 1.0 1.5 2.0 2.5 3.0 ' 0.0 0.5 1.0 1.5 2.0 2.5 3.0 

kAX kAx 



Fig. 3. (Left above) Stability analysis for mixed schemes. The macroscopic variables are set as 
p = 1.0, T = 1.0, ui = 10.0, U2 = 0.0, the constants are set as ci = 10, C2 = 20,770 = 10, At = 
aAx/c2- 

Fig. 4. (Right above) The influence of C2- The macroscopic variables are set as p = 1.0, T = 
1.0, til = 10.0, «2 = 0.0, the other constants are set as ci = 10,r)o = 10, At = 0.45Aa;/c2. 

dissipation. The latter corresponds to the conventional first-order upwind scheme. 
For the given parameters, when the time step is small, both treatments give stable 
simulations; but when the time step becomes large, the first treatment makes the 
simulation stable while the second one does not. 

Figure 4 shows an investigation to the influence of constant C2 on the stability 
of LB simulation, where the value of C2 is altered from 10 to 50, the time step 
At — 0A5Ax/c2, the other constants and macroscopic variables arc unchanged. 
The LB is stable for all tested values of C2. Our experience shows that the value of 
C2 does not influence much the numerical stability if it is not smaller than 2ci , but 
the stable time step becomes smaller for larger the vahie of C2 . Figure 5 shows an 
investigation to the influence of the value of rjQ. The value of 770 is altered from 5 
to 20, C2 = 20, At = 0.3Ax/c2, the other constants and macroscopic variables are 
kept unchanged. We get an indication that it is not difficult to find an appropriate 
value of r]o to get a stable simulation. For cases shown in the figure, only a too small 
value of r]o may result in instability (see the case of r/o = 5) and stability is nearly 
the same when r/o exceeds some critical value (see the cases with 779 = 15 and with 
Vo = 20). 

Since the density p can be normalized to 1, we then investigate only the effects 
of the other two physical quantities, temperature T and fiow velocity u. Figure 6 
shows three cases with different temperatures, T = 1, T = 5 and T = 25. When 
other parameters arc fixed, the numerical stability increase with the increasing 
of the system temperature. This can also be understood that higher temperature 
corresponds to higher sound speed and lower Mach number. 

Figure 7 shows cases with difference flow velocities. The value of ui is altered 
from zero to 15 and U2 = 0. For parameters used in this case, we can flnd that the 
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Fig. 5. (Left above) The influence of rjo. The three physical quanties are set as p = 1.0, T = 
1.0, ui = 10.0, U2 = 0.0, and other constants are set as ci = 10, C2 = 20, r]o = 10, At = 0.3Aa;/c2. 

Fig. 6. (Right above) The influence of temperature T. The other physical quantities are set as 
p = 1.0, «i = 10.0, U2 = 0.0, the constants are set as ci = 10, C2 = 20, rjo = 10, At = 0.45Ax/c2. 




Fig. 7. (Left above) The influence of flow velocity ui. The other physical quantities are set as 
p = 1.0, T = 1.0, U2 = 0.0, the constants are set as ci = 10, C2 = 20, t^q = 10, At = 0.3Ax/c2. 

Fig. 8. (Right above)The x dependence of p, p, u and T. The symbols are simulation results by 
the new LB and lines are analytic solutions. The initial condition is described by Eq. I I25I I. Here 
Ax = 0.01, At = 0.00008, ci = 25,C2 = 50,r?o = 30, terminal time t = 0.36. 



simulation will not be stable if ui is much larger than Ci, even though |a;|maa; is 
only slightly larger than 1 at fcAa; w 0.5. Our experience shows that the value of ci 
can be set nearly equal to the maximum of the flow velocity. 

In summary, constants ci, C2 and 770 influence heavily the stability. In practical 
simulations, ci can be set approximately equal to the maximum of flow velocity; C2 
can be set to be about 2^3 times of the value of ci; r/o can be set an appropriate 
value in between ci and C2. Equation (|24|) is indicative in choosing parameters for 
stable LB simulations of high-speed flow. 
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4. Numerical validations 

Two kinds of benchmarks are used to validate the proposed scheme. The first one 
is the Riemann problem. The second is the problem of shock reflection. 

4.1. Riemann problem 

Here the two-dimensional model is used to solve the one-dimensional Riemann prob- 
lem. The initial macroscopic variables at the two sides are pL, PL and ul, and pR, 
Pn and M/j, respectively. We firstly simulate a Riemann problem with an initial 
condition described by 



where the subscripts "L" and "R" denote the left and right sides of the discontinuity. 
The initial Mach number of left flow is equal to 30.0. The numerical results for 
7 = 1.4 are shown in Fig. 8, where the symbols are simulation results and solid lines 
are analytical solutions. The parameters used in the simulation are ci = 25, C2 = 50, 
j^o = 30, t = 0.36. The size of grid is Ax = Ay = 0.01. Time step Ai = 0.00008. The 
two sets of results have satisfying agreement. In this case no evident "wall-heating" 
phenomenon is observed. As a comparison, we show a result with the general first 
order upwind scheme for the pressure in Fig. 9(a). A abrupt decrease in pressure 
around x = 12 corresponds to the well known wall-heating phenomenon. In order 
to observe the effects of various additional viscosity, we vary the value of Aq from 
ciAx to 5ciAx under the fixed Ai(i = 1,...,8). Figure 9 (b) shows the simulation 
results and the exact one. We can find that the numerical width of shock becomes 
wider and wall-heating problem becomes more pronounced as Aq increases. Results 
in Fig. 9 confirm that (|24p is indicative in choosing the additional viscosity. 

The second example is the propagation of a shock with high ratios of density 
and pressure. The initial macroscopic variables are give by 



The size of grid is Ax = Ay = 2.5 x lO^'^. The numerical results are shown for 7 = 
1.4 in Fig. 10, where the symbols are simulation results and solid lines correspond 
to exact solutions. We find also a good agreement between the two sets of results. 

4.2. Shock reflection 

We will present two gas dynamics simulations. Both are done on rectangular grid. 
The first is to recover a steady regular shock reflection. The second test problem is 
the double Mach reflection of a shock off an oblique surface. This example is used 



PL = 1.4, PR = 1.4, 
PL = 1-0, PR = 1.0, 
UL = 30.0, UR = 0.0, 



(25) 



PL = 1000.0, PR = 1.0, 
PL = 1000.0, PR = 1.0, 
UL = 0.0, UR = 0.0, 



(26) 
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Fig. 9. Comparison of various finite-difference schemes (a) and artificial viscosities (b). The initial 
condition is same as Fig.8. Here Ax = 0.01, At = 0.00008, ci = 25,C2 = 50,r)o = 30, i = 0.36. 
(a) Profile of pressure. The values of Xi(i = 0, ...,8) are the same as those in Fig.8, while /3 = 
0. The symbols correspond to simulation result the line is for analytic solution, (b) Profile of 
temperature. The dashed, dash dotted and dotted lines correspond to Aq = ciAx, 2ciAa;, and 
5ciAx, respectively. The values of \i(i = 1, ...,8) and /3 are the same as in Fig.8. 




in Ref. as a benchmark test for comparing the performance of various difference 
methods on problem involving strong shocks. 

In the first test problem, we have performed a 30° shock reflection for 7 = 1.4. 
The computational domain is a rectangle with length 9 and height 3 (See Fig. 11 (a)). 
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Fig. 11. (See Figll.jpg ){Color online) Regular shock reflection, (a), Sketch map of the steady 
state regular reflection problem, (b), The density contour at time t = 2.5 with Ax = Ay = 0.01, 
At = 1.5 X 10~*, ci = 10, C2 = 20,r;o = 15; Left and up boundary conditions are given by Eg. 1271 1. 
From black to yellow, the value increases. 

Fig. 12. (See Figl2.jpg) (Color online) Double Mach reflection, (a) initial conflguration; (b) the 
density contour at time i = 1.5 with Ax = Ay = 0.01, At = 2.0 X 10~*, ci = 8,C2 = 16,r;o = 10. 
The reflecting wall begins at 20 mesh length from the lower left corner. From black to yellow, the 
value increases. 

This domain is divided into a 900 x 300 rectangular grid with Ax — Ay — 0.01. 
The boundary conditions are composed of a reflecting surface along the bottom 
boundary, supersonic outflow along the right boundary, and Dirichlet conditions on 
the other two sides, given by 

(p, Ui,U2,p)\o,y,t - (1.0, 10.0, 0.0, 1/1.4) 

(p, ui, U2,p)U,i,t = (5.0, 8.0, -3.4641, 20.7143) ^ ^ 

Initially, we set the solution in the entire domain to be that at the left boundary, 
the corresponding Mach number is 10.0. In Fig. 11(b) we show a contour plot of the 
density. The clear shock reflection on the wall agrees well with the exact solution. 

The second test problem is an unsteady shock reflection. A planar shock is 
incident on an oblique surface with the surface at a 30° angle to the direction of 
propagation of the shock (Fig. 12(a)). The fluid in front of the shock has zero velocity, 
and the shock Mach number is 10.0. In Fig. 12(b) we show the result of density 
contour, where the double Mach reflection phenomenon is successfully recovered. 

5. Conclusions and discussions 

The lattice Boltzmann simulation to high-speed compressible flows is revisited by 
proposing an improved LB model. The new LB model is composed of the origi- 
nal discrete-velocity-model by Kataoka and Tsutahara and an appropriate finite- 
difference scheme to the convection term. An additional dissipation term is intro- 
duced to improve the numerical stability. The adding of the dissipation term should 
survive the dilemma of stability versus accuracy. In other words, the dissipation 
should be minimal but make the evolution satisfy the von Neumann stability con- 
dition. The effects of polynomial equilibria-'^^ are taken into account (via the first 
term of Eq. (|23p ) in such an approach. Due to the complexity the analysis resorts to 
the software, Mathematica-5, and only some typical results are shown by figures. 

Benchmark tests are used to validate the proposed scheme and reference values 
of model parameters are suggested. Typical Riemann problems with high-Mach- 
number (30 or higher) and high ratios (1000 : 1) of pressure and density show good 
accuracy and stability of the new scheme, even though they are generally difficult 
to resolve by traditional computational fluid dynamics. Regular and Mach shock 
reflection problems are successfully recovered, which shows also the potential appli- 
cation of lattice Boltzmann model to fluid systems where non-equilibrium processes 
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are intrinsic and pronounced. The new LB model may be used to investigate some 
long-standing problem, such as the transition between regular and shock reflections. 
At the moment, wc arc still not able to present a complete description on the most 
appropriate additional dissipation term, but the idea presented in the paper can be 
easily used to get some practically useful solutions for stability enhancement. We 
plan to better clarify the physical dissipation and artiflcial ones in the future. 
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